Efficient DMFT-simulation of the Holstein-Hubbard Model 
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We present a method for solving impurity models with electron-phonon coupling, which treats the 
phonons efficiently and without approximations. The algorithm is applied to the Holstein-Hubbard 
model in the dynamical mean field approximation, where it allows access to strong interactions, 
very low temperatures and arbitrary fillings. We show that a renormalized Migdal-Eliashberg theory 
provides a reasonlable description of the phonon contribution to the electronic self energy in strongly 
doped systems, but fails if the quasiparticle energy becomes of order of the phonon frequency. 
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Electron-lattice interactions play a fundamental role 
in the physics of metals. In conventional ("weakly cor- 
related") materials such as Al or MgB2, phonons are 
the dominant source of inelastic scattering and provide 
the pair binding which leads to superconductivity. In 
these materials, physical understanding is greatly aided 
by the Migdal-Eliashberg theory, which justifies the use 
of density functional band theory methods to estimate 
the electron-phonon coupling constants and perturbative 
methods to estimate electron mass renormalizations, life- 
times and pairing strengths. A key result of the Migdal- 
Eliashberg theory is that the effect of the electrons on 
the phonons is a renormalization of the phonon oscilla- 
tor frequency ojq, while the effect of the phonons on the 
electrons is an increase in the electron effective mass. 
The mass increase "turns on" at frequencies below the 
renormalized phonon frequency. 

In unconventional ( "strongly correlated" ) materials 
such as the high temperature superconductors [l|, the 
fullerenes @, El], or the colossal magnetoresistance rare- 
earth manganites 0, H| the situation is less clear. Many 
experiments suggest that electron-phonon effects are 
important. For example, in high T c superconductors 
changes attributed to the electron-phonon coupling are 
observed in the electronic dispersion [l[ at energies of 
the order of typical optical phonon frequencies. How- 
ever, there is as yet no clear theoretical basis for inter- 
preting these or related data in strongly correlated mate- 
rials. While there has been extensive and important work 
on static properties of models involving electron-electron 
and electron-phonon coupling @, 0] , on Fermi- liquid and 
effective field theory based approaches [1, Q , on the Hol- 
stein model [l(j and on the ( bi-)p olaron problem (one or 
two interacting electrons) [ll|, [l2| less is known about the 
dynamical consequences of the electron-phonon interac- 
tion in strongly correlated materials. 

A way forward is provided by dynamical mean field 
theory (DMFT), a powerful, non-perturbative tool to 
study the properties of strongly correlated systems, 
which has provided considerable insight into the cor- 
relation induced metal- insulator (Mott) transition [3]. 
Progress in using dynamical mean field methods to 
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bl . la . 1^ 1 has been hampered by the mismatch in en- 
ergy scales between the electronic and lattice parts of the 
problem and by the difficulty of providing a quantitative 
treatment of the low temperature properties of strongly 
correlated materials even in the absence of electron- 
phonon coupling. In this paper we introduce a new 
method which resolves these problems and allows the 
DMFT simulation of important classes of models at es- 
sentially the same computational expense as the corre- 
sponding models without coupling to phonons. 
We consider the Holstein-Hubbard Hamiltonian 
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where U denotes the on-site repulsion, /i the chemical 
potential of the electrons with creation operators cj. and 
density operators n a , the creation operator for Ein- 
stein phonons of frequency liq, and the electron-phonon 
coupling is A. 



The single-site DMFT approximation [13[ reduces the 
problem to the self-consistent solution of a quantum 
impurity model specified by the Hamiltonian Hqi = 
H\ oc + i/hyb + -ffbath- Here, the local term is 



H\ oc = -fJ>(n-\ + "-1) + Unyni 

+A(n T + n l - l)(6 t + b) + uj^b, 



(2) 



and the impurity-bath mixing and bath Hamiltonians 
are H hyh = J2 P ,a( v p,<? c t a P,o- + Vp,o- c <? a i,o-) and # b ath = 



Tip a e p a p o a p-o- The parameters V v , a and e p are deter- 
mined by a self-consistency equation. 

In the absence of electron-phonon coupling, H\ oc has a 
small Hilbert space and one energy scale (U). Adding an 
electron-phonon coupling introduces a new energy scale, 
the phonon frequency uj , and requires keeping track of 
the bosonic sector of the Hilbert space (with an infinite 
number of states). Previous approaches to the problem 
have involved either treating the bosons semiclassically 
or truncating the boson Hilbert space, retaining 
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only a finite number of boson states 15|, [l6| . The semi- 
classical approach cannot account for quantal phonon ef- 
fects such as electronic mass renormalization or supercon- 
ductivity, while treating even a truncated boson Hilbcrt 
space directly is computationally expensive. 

In Refs. [3, [l9[ we have shown that the stochastic 
sampling of a diagrammatic expansion of the partition 
function in powers of the impurity-bath hybridization 
term i/hyb leads to a highly efficient impurity solver 
for purely electronic models. After tracing out the bath 
states a p>a , the weight of a Monte Carlo configuration 
corresponding to a perturbation order n (n creation op- 
erators cJ.(t ct ) and n annihilation operators c ct (t£)) can 
be expressed as 

«>({0ifa)}) = Tr c (T T e--fo dTH ^0 2n (r 2n ) . . . 

...0 2 (t 2 )Oi(ti)) dr 1 ...dT 2 „n(detM- 1 ) S(T , (3) 
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FIG. 1: Phase diagram for uj = 0.2t, fit = 50,400 and 
U/t = 6 in the plane of chemical potential jj, and phonon 
coupling strength A. Half-filling corresponds to /i = U /2. In 
the metallic phase, the filling increases with increasing fj,, A. 



with Oi{ji) the (time ordered) creation and annihilation 
operators for spin up or down electrons on the impurity 
site. The matrix elements (M~ 1 )ij — F^^t^ i — T a.j) are 
determined by the V Pl(T and e p through the hybridization 



functions Fa-(-iu n ) = J2 P 



[19|. The sign s a is 1 



if the cr-spin operator with the lowest time argument is 
a creation operator and —1 otherwise. 

To evaluate (•••)$, we use a Lang-Firsov [2l[ trans- 
formation. Defining operators X = (tf + 6)/V2 and 
P = (b^ — b)/iy/2, the unitary transformation specified 
by e iPXo shifts IbyX = {V^X/too)^ +n l - 1) so that 



ffioc = e iPXo H loc e 



-iPXo 



= -f l (h ]+ n i ) + Un^n i + ^{X 2 + P 2 ) (4) 

has no explicit electron-phonon coupling. H\ oc is of 
the Hubbard form but with modified chemical poten- 
tial and interaction strength: jl = /j, — X 2 /luq, U — 
U — 2X 2 /llIq- Also, the electron creation and annihilation 
operators are transformed according to c£ = b ^c£ 

and c a = e ~^e CT . The phonon contribution Wf, to 
the weight © is w 6 ({O l (r l )}) = ( e s ^^(^„) . .. e «iA( n )^ 
with < n < . . . < Ti n < S{ = 1 (-1) if the 
n th operator is a creation (annihilation) operator and 
A(t) = ^(e^tf - e^° T b). This expectation value, 
taken in the thermal state of free bosons, evaluates to 



Wb{{Oi{Ti)}) = exp 



A 2 K 2 
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1) 



n-r,) } ) 



,(5) 



2n>i>j>l 

and the weight ([3]) becomes the product 

w({0<(n)}) = «7i({O i (ri)})tSHubb«d({O < (r0}), (6) 



where WHubbard is the weight of a corresponding configu- 
ration in the Hubbard impurity model (without phonons, 
but modified parameters U and jX). The Holstein- 
Hubbard model can thus be simulated without trunca- 
tion at an expense comparable to the Hubbard model. 
Superconducting phases can be simulated in the same 
way, using the formalism outlined in Ref. [l3j ]: the only 
difference is that the determinant in Eq. ([3]) can no longer 
be factorized into spin components. 

We have applied our method to the Holstein-Hubbard 
model with a semi-circular density of states of bandwidth 
At and uio/t = 0.2. For reasons of space, we restrict 
our attention to non-superconducting phases. Figure [1] 
shows the phase diagram at U/t = 6 in the space of 
chemical potential and phonon coupling. Our results at 
half-filling are in good agreement with those computed 
by NRG [15|. If the chemical potential is increased at 
fixed electron-phonon coupling, a transition to a metal 
occurs. The shift in critical /i with temperature arises 
from the entropic stabilization of the insulating phase 
and is roughly A- independent. The transition is first 
order at T > (with a jump in density), but appar- 
ently becomes continuous as T — > 0. The critical chem- 
ical potential decreases as the electron-phonon coupling 
is increased. This behavior is physically expected: in- 
creasing the phonon coupling reduces the effective inter- 
action U and so the magnitude of the insulating gap. We 
have confirmed this by analysis of the Green function. 
Our finding appears to contradict Ref. [n| which stated 
that the electron-phonon coupling stabilizes the insulat- 
ing state. The apparent difference arises from a choice of 
convention: the authors of Ref. fl6[ couple the phonons 
to the total density n, which implies a A-dependent shift 
in chemical potential which was interpreted as a physical 
stabilization. We couple the phonons to the difference in 
density from the half filled value, which eliminates this 
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(physically irrelevant) shift at half filling. 

We now turn to the frequency dependence of the elec- 
tron self energy. In weakly correlated materials the 
conduction-band electrons renormalize the phonon prop- 
agator D from its bare value Dq according to the RPA 
(ladder) result D^ 1 = D^ 1 — X 2 Xo with \o the bare 
(no phonons) electronic density-density correlator. The 
Migdal-Eliashberg approximation also implies that the 
local density-density correlator \ is given by 

X (0) = xo(0)/(l - A 2 A)(tt)xo(N)). (7) 

The electron self energy is given in frequency space by 
the convolution £ = X 2 Gq*D. It is convenient to present 
results as the Matsubara axis mass renormalization fac- 
tor r(iu) n ) = —^sm^(iuj n )/uj n . According to Migdal- 
Eliashberg theory, the effect of the electron-phonon cou- 
pling is to add to the purely electronic r{iuj n ) a term 
approximately of the form 

r ME {iUn) = (2A 2 A r o/w n )arctan(w n /wo cn ) (8) 

with No the Fermi surface density of states (in our case 
l/2irt). The zero frequency limit gives the electron- 
phonon contribution to the Fermi surface mass enhance- 
ment m* ME /m — 2 ^,-£f° and r ME drops to about half of 
its zero frequency value around 2o;o cn . 

The upper panel of Fig. [2] compares the density-density 
correlation function computed from Eq. ([7]) to x(r) mea- 
sured in our QMC simulations, for U = and several 
values of A. Good agreement is seen; the small shift of A 
required to match the A = 0.3 QMC data (A c ff ~ 0.275) 
is a beyond-Migdal effect arising from the small but non- 
vanishing effect of the electron self energy on xo- x( T ) m " 
creases with A because density fluctuations can transform 
into phonons. The intermediate-time exponential decay 
visible in the figure gives the renormalized phonon fre- 
quency ujQ Cn , while in the regimes where there is no clear 
evidence of exponential decay the phonon frequency is 
not significantly renormalized. For U > the computed 
x(t), shown in the lower panel of Fig. [21 also exhibits an 
exponential decay, from which we estimate u>Q en . 

Figure [3J presents the phonon-induced change 
in the mass renormalization function Ar = 
— Sm (S(A, 8, u> n ) — E(0, 8, oj n )) /u> n . The upper panel 
shows results for U = it, less than the Mott critical value 
U C 2 ~ 5.8t. At half filling (curves traced out by symbols) 
we find, in agreement with Ref. [H|], that the phonon 
contribution to the mass enhancement has the "wrong" 
sign, except for couplings very close to the bipolaronic 
instability: the dominant effect of the phonons is the 
reduction of U. A new result is that (again except for the 
largest A) there is no obvious feature at the renormalized 
phonon frequency. The blue lines present the doping 
dependence at fixed A = OAt. We see that as the doping 
is increased, a result compatible with Migdal-Eliashberg 
theory is eventually recovered. The lower panel shows 




FIG. 2: Imaginary time dependence of the density- density 
correlation function plotted on a log-scale to highlight the 
intermediate time exponential decay. If no exponential regime 
is visible, the phonon frequency is not renormalized; otherwise 
u>o cn can be extracted from the slope. Upper panel: U — 
and n = 1. Dashed lines: xij) from Eq. for indicated 
phonon couplings. Solid lines with 1-a error bars show the 
QMC results for fit = 400 and X/t = 0.1, 0.2, 0.3 (from 
bottom to top). Lower panel: correlation functions for the 
interacting model at half filling, U/t = 4 and X/t = 0.6, 0.65 
and in the doped Mott insulator at U/t — 6, X/t — 0.4, 0.6 
(dopings per spin 8 = 0.043, 0.138). 



the behavior for U = 6t, greater than the Mott critical 
value. Here we find again that for moderate electron 
phonon couplings and relatively highly doped samples, 
a renormalized Migdal-Eliashberg theory (with perhaps 
a slightly reduced effective electron phonon coupling 
A e ff) provides a reasonably consistent description of the 
data, except that Ar drops more rapidly for u> n ^> luq- 
However, when the quasiparticle energy Ep /(l — dY>/duo) 
becomes of order of the phonon frequency (as happens 
at 8 = 0.04) or the coupling approaches the bipolaronic 
instability (A = 0.6t) the Migdal-Eliashberg description 
breaks down. For large A, rME{iu n ) drops too rapidly, 
while near the doping driven Mott transition, the shape 
is too broad. Nevertheless, unlike in the half-filled metal 
shown in the upper panel of Fig. [31 in the weakly doped 
Mott insulator, the upturn in the mass renormalization 
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FIG. 3: (Color online) Phonon contribution to the elec- 
tronic self energy at j3t = 400. Black diamonds: renormalized 
phonon frequency ujQ Cn . Upper panel: U/t — 4. Symbols: half 
filling, A values indicated. Blue lines: \/t — 0.4 and doping 
per spin 8 = 0.019, 0.067 and 0.182 (from bottom to top). 
Black line: Migdal-Eliashberg result. Lower panel: U = 6t 
(above the critical value for Mott insulating behavior). Sym- 
bols: A and 8 as shown. Solid lines: Migdal-Eliashberg calcu- 
lations, Eq. (JHJ, with renormalized electron phonon couplings 
A ren and phonon frequencies cjQ Cn as indicated. 



still coincides with w^ 8 ". Thus, although a Migdal- 
Eliashberg analysis of self energies in high-T c materials 
is problematic, the association of kinks in dispersions 
with phonon features may be robust. 

To summarize: we have introduced a method which 
enables simulations of impurity models with electron- 
phonon coupling at essentially the same computational 
cost as the corresponding models without phonons. We 
used the method to compute the phase diagram for the 
doping driven Mott transition in the Holstcin-Hubbard 
model and the frequency dependence of the self energy. 
Previous work by one of us and by other groups [1, [lj| 
had argued on the basis of Fermi liquid theory that the ef- 
fect of phonons on correlated systems should be similar to 
that in uncorrelated systems, but with renormalized cou- 
plings. We found that this expectation is simply wrong 
at strong interactions in a half-filled model, but applies 
for large enough fillings and in the doped Mott insulating 



state as long as the quasi-particle energy is larger than 
the phonon frequency. 

The approach presented here generalizes to other mod- 
els which can be decoupled by Lang-Firsov transforma- 
tions. In the Holstein-Hubbard model the interplay be- 
tween superconductivity and strong correlations is an im- 
portant open problem. Further generalizations are possi- 
ble, for example to Jahn- Teller couplings in multiorbital 
models - although if pair-hopping terms are important 
one must perform a double-expansion in these and in 
the hybridization. Application of the method introduced 
here to models of manganites and to Cqq (if the elec- 
tronic spectrum is truncated to the 5 levels nearest to 
the chemical potential) seems entirely feasible. 

The calculations have been performed on the Hreidar 
beowulf cluster at ETH Zurich, using the ALPS-library 
[23I ]. We thank M. Troyer for the generous allocation 
of computer time, K. Ziegler for helpful discussions, and 
NSF-DMR-040135 for support. 
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